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Abstract 

The ground-state phase diagram of the asymmetric Hubbard model is stud- 
ied in one and two dimensions by a well-controlled numerical method. The 
method allows to calculate directly the probabilities of particular phases in 
the approximate ground-state and thus to specify the stability domains corre- 
sponding to phases with the highest probabilities. Depending on the electron 
filling n and the magnitude of the asymmetry tf/td between the hopping inte- 
grals of / and d electrons two different scenarios in formation of ground states 
are observed. At low electron fillings (n < 1/3) the ground states are always 
phase segregated in the limit of strong asymmetry (td ^> £/)• With decreasing 
asymmetry the system undergoes a transition to the phase separated state 
and then to the homogeneous state. For electron fillings n > 1/3 and weak 
Coulomb interactions the ground state is homogeneous for all values of asym- 
metry, while for intermediate and strong interactions the system exhibits the 
same sequence of phase transitions as for n small. Moreover, it is shown that 
the segregated phase is significantly stabilized with increasing electron filling, 
while the separated phases disappear gradually from the ground-state phase 
diagrams. 

PACS nrs.:75.10.Lp, 71.27.+a, 71.28.+d, 71.30.+h 



1 



1 Introduction 



The asymmetric Hubbard model is one of the simplest models for a description of 
correlated fermions on the lattice. It has been used in the literature to study a 
great variety of many-body effects in rare-earth and transition-metal compounds, of 
which quantum phase transitions, mixed-valence phenomena, charge- density waves, 
and electronic ferroelectricity are the most common examples [1, 2, 3, 4, 5]. In the 
last few years the asymmetric Hubbard model was also used for a description of 
ground-state properties of fermionic particles on the optical lattice [6, 7]. The model 
consists of two species of electrons: heavy / electrons and light d electrons. The 
Hamiltonian of the model is 

H = -t d £ dtd s - t f £ ftfj + U £ ft fifth (l) 

<ij> <ij> i 

where /+ (/j) and df (dj) is the creation (annihilation) operator of heavy and light 
electron at lattice site i. 

The first two terms of (1) are the kinetic energies corresponding to quantum- 
mechanical hopping of d and / electrons between the nearest neighbor sites i and 
j with hopping probabilities td and tf, respectively. The third term represents the 
on-site Coulomb interaction between the d electrons with density rid = jY^i dfdi 
and the / electrons with density rif = jJ2iftfi, where L is the number of lattice 
sites. The model is called "asymmetric" because the hopping integrals for d and / 
electrons may be different. Usually, the hopping integral of the d electrons is taken 
to be the unit of energy (td = 1) and the /-electron hopping integral is considered 
in the limit ^< 1. This is a reason why the d electrons are called light and the / 
electrons heavy. The Hamiltonian (1) reduces to the spinless Falicov-Kimball model 
for tf — and to the usual one-band Hubbard model for tf = 1. Thus the asymmetric 
Hubbard model can we viewed as a generalized Falicov-Kimball model and also as a 
generalized one-band Hubbard model. 

The first systematic study of ground-state properties of the asymmetric Hubbard 
model has been performed by Lyzwa et al. using various analytical and numerical 
techniques [1, 8, 9, 10]. In the first paper from this series the authors studied the 
ground-state properties of the one dimensional asymmetric Hubbard model by own 
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approximate method that allowed to treat larger clusters than accessible by exact 
diagonalization technique. Their method was based on the sequence of two steps. 
First, they found the lowest-energy state for every permissible /-electron configura- 
tion, similarly as in the pure spinless Falicov-Kimball model (tf — 0). To do this 
the matrices of rank ^ L _^' d y Nd \ had to be diagonalized. Second, they took the states 
thus found as a basis of a new matrix of rank ^ L _^ f y Nf \ that was subsequently di- 
agonalized. The lowest-energy eigenstate of this matrix was then used to construct 
the approximate ground state. The main result obtained by the application of this 
method to the asymmetric Hubbard model was that the motion of the heavy elec- 
trons is strongly influenced by the light ones, while the light electrons are almost 
unaffected by the presence of the heavy ones. The subsequent study [8] of the asym- 
metric Hubbard model on the one- dimensional ring with two / and two d electrons 
showed that an effective attraction between two /-electrons can be produced by cor- 
relation effects for a certain set of the model parameters. The same result, an effective 
attraction between two heavy electrons mediated by two light electrons leading to the 
phase segregation, was confirmed also in two dimensions [10]. These studies showed 
that the phase segregation is own not only to the Falicov-Kimball model [11], but 
persists also at finite tf. For strong interactions this result was proven rigorously 
by Ueltschi [12]. The boundary of the phase segregation/separation region in the 
U — tf plane has been calculated very recently by two different methods. To identify 
the transition boundary Gu at al. [7] used the quantum entanglement between a 
local part and the rest of the system and the structure factor of charge- density wave 
(CDW) for heavy electrons. Away from half-filling, they found that the domain of 
phase separation is not confined only to small tf, but persists up to relative large 
values of tf (e.g., tf ~ 0.2 for intermediate Coulomb interactions U ~ 5) and with 
increasing U is further stabilized. The same result has been obtained also by Wang 
et al. using the bosonization method [5]. 

In the current paper we study the ground-state phase diagram of the asymmetric 
Hubbard model by an improved Lyzwa's scheme [1] discussed above. The advantage 
of this method is that it can treat much larger clusters than accessible by exact diag- 
onalization technique and its applicability, unlike the DMRG [7] or bosonization [5] 
method, is not confined only to the one-dimensional case. Moreover, the method 
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allows to calculate directly the probabilities of particular / electron configurations 
and thus to specify the stability domains corresponding to distributions with the 
highest probabilities. 

2 Method 

Before discussing our approach, it is useful and instructive to summarize main steps of 
the numerical algorithm leading to the exact solution of the spinless Falicov-Kimball 
model on finite clusters. The asymmetric Hubbard model (1) reduces to the spinless 
Falicov-Kimball model for tf — 0. Its Hamiltonian reads 

H FKM = -U E d+d, + ftfidtdi. (2) 

<ij> i 

Since in this version, without /-electron hopping, the /-electron occupation num- 
ber /j + /j of each site % commutes with the Hamiltonian (2), the /-electron occupation 
number is a good quantum number and can be replaced by classical variables Wi = 1 
or 0, according to whether or not the site i is occupied by the localized / electron. 
Then the Hamiltonian (2) can be written as 

^irKM = ^hijdfdj, (3) 

ij 

where hij(w) = —td, if i and j are the nearest neighbor; hij(w) = Uwi, if i = j and 
zero otherwise. 

Thus for a given /-electron configuration w = {wi, . . . wl} defined on the one, 
two, or three-dimensional lattice, the Hamiltonian (3) is the second- quantized version 
of the single-particle Hamiltonian and can be directly diagonalized by the following 
canonical transformation 

dtM = T,v^di, (4) 

i 

where is the unitary matrix that diagonalizes h(w). The ground state of H-tkm 
is then constructed as follows 

l<) = EKHlo) (5) 

a=l 
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and the corresponding ground-state energy is simply given by 

N d 

E o(w) = E £ «H ' ( 6 ) 

a=l 

where e a {e\ < 62 < ■ ■ ■ < £l) are eigenvalues of the single particle matrix h(w). 

Generalizing this procedure our approach to the full Hamiltonian of the asym- 
metric Hubbard model (1) can be formulated in the following two points: First, we 
construct the reduced basis \ipk) of H by making the Ansatz 

m = \4M), (7) 

where \tp[) is the complete set of eigenstates of the /-electron subsystem (k = 
1,2, ... , (j,_jvy)jjv- ; ) an d iV'fc) is the ground state corresponding to \ip[) (note that 
ft fi\i , t) =w ^\'4 , l))- Second, the reduced basis \ipk) is used to calculate the matrix 
elements H nm = (xjj n \H\xjj m ) of the full Hamiltonian (1). This matrix is then diago- 
nalized and its lowest energy eigenvalue E\ yields the upper bound for the ground 
state energy of H. The corresponding eigenvector 

hfoHE^.ihM > (8) 

n 

(where U nm is the unitary matrix that diagonalizes H nm ) is the approximate ground 
state and can be used directly to calculate the expectation value of any operator 
(A) = {^q\A\^q) . Using expressions (4-8) one can show easily that any ground- 
state expectation value can be written directly in terms of the unitary matrices 
that diagonalize the single-particle matrices h{w) corresponding to all possible 
distributions of /-electrons. 

To establish connection between our approach and Lyzwa's one let us calculate 
explicitly the matrix elements H nm of the total Hamiltonian of the asymmetric Hub- 
bard model (1). Making use of the fact that the new creation and annihilation 
(d a ) operators obey the following commutation relation 

d a (n)d^(m) + d^(m)d a (n) = k a p(n,m) , (9) 

where the matrix elements k a p are given by 

Mn>™)=EV& B) V& m) , (10) 
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the straightforward calculations lead to the following expressions for the diagonal 
and off-diagonal matrix elements of H 

H nn = E (n) (11) 

and 

H nm = -t f (^\^ftfM) = -(-m f det(k(n,m)) , (12) 

(ij) 

where s is the number of permutations that should be done to transform ip^ on ip^ 
by fffj. As discussed bellow the Lyzwa's approach can be recovered directly from 
Eq. 12 by putting det(k) = 1, however generally det(k) ^ 1. 

To test our method we have first calculated the ground-state energy of the asym- 
metric Hubbard model for various model parameters (tf, U) on small one-dimensional 
clusters where exact results are also accessible. In Table 1 and Table 2 we present 
the exact and approximate ground-state energies obtained for the finite clusters of 
L = 6 and 10 sites, two representative values of tf (tf— 0.1, 1) and three represen- 
tative values of U ({7=0.1 , 1, 10). For a comparison we have included into Table 1 
also results obtained by Lyzwa et al. [1] and we have verified numerically that these 
results can be recovered exactly by our method simply putting det(fc) = 1 in the 
Eq. 12. Comparing these results one can see that both approaches work very well 
in the weak-coupling limit (U < 1), while in the opposite limit (U >> 1) our method 
yields a considerable better estimation of the ground-states energy than the Lyzwa's 
one. Fig. 1 demonstrates that this trend holds also for smaller values of tf. More- 
over, as one can expect intuitively, the accordance between our results and the exact 
ones considerably improves with decreasing tf and for a sufficiently strong asym- 
metry (tf ~ 0.2) a nice accordance of results is observed over the whole interval of 
Coulomb interactions (see also Table 2). 

To verify the ability of our method to describe the main characteristics of the 
exact ground state we have also calculated the /-electron pair correlation function 
L(x) defined by [1] 

/-(•'•) {yiififj-fj-J,-,)- (13) 

j 

The results of numerical computations are summarized in Fig. 2. It is seen that a 
nice accordance with exact behaviors is obtained in the weak (U = 1) as well as 
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strong (U = 10) coupling limit for both small (tf < 0.2) and intermediate (tf ~ 0.4) 
values of /-electron hopping integrals. 

3 Phase diagrams 
3.1 One-dimensional case 

One of the greatest advantages of the method discussed above is that it allows to 
calculate directly the probability of any /-electron configuration in the approximate 
ground state. We used this fact to construct the phase diagrams of the asymmetric 
Hubbard model in the tf — U plane for various sizes of clusters and electron fillings. 
The phases presented in the phase diagrams are those corresponding to /-electron 
distributions with the highest probabilities for given values of tf and U . 

The typical examples of ground-state phase diagrams of the asymmetric Hubbard 
model are displayed in Fig. 3 for L = 24 and four representative values of /-electron 
filling rif (note that rid = rif = n/2). A general feature which can be noticed in these 
pictures is that the basic structure of the phase diagrams is formed by only two main 
types of the /-electron configurations, and namely, the most homogeneous distribu- 
tions MHD (the / electrons are distributed homogeneously over the whole lattice) 
and the phase separated configurations (the / electrons occupy only one part of the 
lattice while the remaining one is empty). As discussed below, between the phase sep- 
arated configurations the special role play the phase segregated configurations (all / 
electrons clump together) and therefore they are considered here as the independent 
group. In the language of an effective interaction between the / electrons, the most 
homogeneous configurations correspond to an effective repulsion and the phase sep- 
arated/segregated configurations to an effective attraction between the / electrons. 
From this point of view, it is very interesting that in the pure electronic system 
(with only the on-site Coulomb repulsion between the light and heavy electrons) an 
effective attraction between the / electrons is produced, even for tf away from the 
Falicov-Kimball limit tf = 0. Indeed, our results show that the phase boundary t C f(U) 
between the most homogeneous and phase-separated/segregated regions increases 
rapidly with increasing U and reaches the intermediate values t C f(U) ~ 0.25 already 
for intermediate Coulomb interactions. In the weak coupling and low density limit, 
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the phase boundary scales like tf(U) ~ U 2 , while at higher electron fillings (rif = 1/4 
and n/ = 1/3) the phase-separated/segregated distributions are stabilized above 
some critical value of Coulomb interaction U (U ~0.5forn/ = l/4 and U ~ 2.2 for 
n/ = 1/3). In these limiting cases our results reproduce the analytical and numer- 
ical results obtained recently by bosonization [5] and exact-diagonalization/DMRG 
method [7]. Comparing our results with the exact-diagonalization/DMRG results [7] 
(rif = 1/4) one can see that these results agree very well in spite of the fact that 
fully different approaches have been used to identify the phase separated region. 

The advantage of our method is that it also allows us to identify the internal 
structure of the phase diagrams and thereby to study how this structure changes 
by varying the model parameters. Fig. 3 shows that the phase diagrams of the 
asymmetric Hubbard model (strictly said the phase-separated domains) have a rich 
internal structure that exhibits some general trends. First, the phase separated region 
starts with the phase segregated distribution. Small exceptions are found only for 
rif = 1/4 and rif = 1/3, where also some other phases are observed for tf — > 0, but 
their stability regions are very limited. Second, the segregated cluster of length Nf 
splits into two or more smaller clusters with increasing tf. Third, the segregated 
configuration is stabilized with increasing electron filling, while the separated phases 
disappear from the phase diagrams. 

Although the cluster used in our numerical calculations is relatively large (L = 24) 
to exclude completely the finite size effects the same calculations have been performed 
on several different clusters for each selected value of electron filling. We have found 
that the fundamental characteristics of the phase diagrams discussed above, and 
namely, the phase boundary between the phase-separated and most homogeneous 
phase, the phase boundary of segregated phase, the critical value of Coulomb inter- 
actions at which the phase separation starts for large /-electron fillings are almost 
independent of L. This is clearly demonstrated in Fig. 4 where the boundary of phase 
separation is plotted for two different values of nj (rif = 1/4 and rif = 1/3). This 
analysis indicates that our one-dimensional results can be extrapolated satisfactorily 
to the thermodynamic limit (L — > oo). 



8 



3.2 Two-dimensional case 



Since the phenomenon of phase separation is one of the most interesting problems in 
the condensed matter physics we extend our calculations also on the two-dimensional 
case. As discussed above such an extension is possible due to the fact that instead 
the full Hilbert space we work only with the reduced basis and corresponding reduced 
matrices of rank ^ L _^ f y Nf \ ■ This allows us to study the two-dimensional clusters up 
to 6 x 6 sites that are far away beyond the reach of present day computers within 
the exact diagonalization calculations. The results of our numerical calculations 
obtained for two different values of electron fillings are summarized in Fig. 5 in 
the form of tf — U phase diagrams together with the complete list of /-electron 
configurations with the highest probabilities. It is seen that all main features of 
the one-dimensional phase diagrams hold in two dimensions, too. For small values 
of /-electron hopping the system is phase segregated/separated, while increasing 
tf stabilizes the homogeneous distribution of / electrons. In accordance with the 
one-dimensional case we have found that the phase boundary t C f{U) between the 
homogeneous and phase separated region scales like U 2 for weak interactions, and 
that the region of phase segregation/separation increases rapidly with increasing rif. 

In summary, we have presented an improved numerical scheme for calculating 
ground-state properties of the asymmetric Hubbard model. The advantage of this 
method is that it can treat much larger clusters than accessible by exact diagonal- 
ization technique and its applicability, unlike the DMRG or bosonization method, 
is not confined only to the one-dimensional case. Moreover, the method allows to 
calculate directly the probabilities of particular / electron configurations and thus to 
specify the stability domains corresponding to distributions with the highest prob- 
abilities. We have used this fact to construct the ground-state phase diagrams of 
the asymmetric Hubbard model in one and two dimensions for wide range of model 
parameters. We have found that at low electron fillings (n < 1/3) the ground states 
are always phase segregated for a strong asymmetry between the hopping integral 
of d and / electrons (t^ ^> tf). With decreasing asymmetry the system undergoes 
a transition to the phase separated state and then to the homogeneous state. For 
electron fillings n > 1/3 and weak Coulomb interactions the ground state (in one 
dimension) is homogeneous for all values of asymmetry, while for intermediate and 
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strong interactions the system exhibits the same sequence of phase transitions as for 
n small. 
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Table 1: The ground-state energy of the one- dimensional asymmetric Hubbard model 
calculated for three different values of Coulomb interaction U on finite clusters of 
L = 6 and L = 10 sites at tf — 1 and iif = = 1/2. Different columns correspond to 
exact results (Exact), Lyzwa's approach (Approx. 1) and our approach (Approx. 11). 





L=6 


L=10 


u 


Exact 


Approx. 1 


Approx. II 


Exact 


Approx. I 


Approx. II 


0.1 


-1.30850 


-1.30868 


-1.30830 


-1.26960 


-1.26978 


-1.26934 


1.0 


-1.10019 


-1.11770 


-1.08064 


-1.06144 


-1.08013 


-1.03697 


10 


-0.27739 


-0.73823 


-0.20421 


-0.27037 


-0.71753 


-0.19908 



Table 2: The ground-state energy of the one- dimensional asymmetric Hubbard model 
calculated for three different values of Coulomb interaction U on finite clusters of 
L — 6, 10 and L = 14 sites at tf = 0.1 and rif — — 1/2. Different columns 
correspond to exact results (Exact) and our approach (Approx. II). 



u 


L 


,=6 


L 


= 10 


L 


=14 


Exact 


Approx. II 


Exact 


Approx. II 


Exact 


Approx. II 


0.1 


-0.70864 


-0.70864 


-0.68725 


-0.68725 


-0.68152 


-0.68151 


1.0 


-0.51544 


-0.51507 


-0.49546 


-0.49495 


-0.49047 


-0.48986 


10 


-0.10056 


-0.09961 


-0.10015 


L -0.09925 


-0.10018 


-0.09923 
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Figure Captions 

Fig. 1. The ground-state energy of the asymmetric Hubbard model as a function 
of U calculated for four different values oft/ and L = 10. Different lines correspond to 
exact results (solid line), our approach (dashed line) and Lyzwa's approach (dashed- 
dotted line). The half-filled band case (rif = n d = 1/2). 

Fig. 2. The /-electron pair correlation function L(x) of the asymmetric Hubbard 
model calculated for two different values of U and tf at L — 10. Different lines 
correspond to exact results (solid line), our approach (dashed line) and Lyzwa's 
approach (dashed-dotted line). The half-filled band case (rif = n d = 1/2). 

Fig. 3. The ground-state phase diagram of the one-dimensional asymmetric Hub- 
bard model calculated for several /-electron densities on finite cluster of L = 24 
sites. Depicted domains represent the stability regions of the /-electron configura- 
tions with the highest probabilities in the approximate ground state. Two small 
domains w a and Wb (for rif = 1/3) consists of several subdomains: w a = [110000]4, 
w\ = [110011000000] 2 ; w\ = [1100] 4 00000000, w 2 b = 111001100111000000000000, 
wl = [111100] 2 000000000000, where the lower index denotes the number of repeti- 
tions of the block [...]. 

Fig. 4. The phase boundary between the homogeneous (MHD) and phase sepa- 
rated (PS) domain as a function of U . Comparison of our numerical results obtained 
on different finite clusters for rif = 1/4 and rif = 1/3. 

Fig. 5. The ground-state phase diagram of the two-dimensional asymmetric Hub- 
bard model calculated for two /-electron densities on finite cluster of L = 36 sites 
with the complete list of /-electron configurations with the highest probabilities. 
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t=0.1 




-0.1 
-0.2 
-0.3 
-0.4 
-0.5 
-0.6 
-0.7 
-0.8 























Exact 






Lyzwa's approach 






Our approach 



4 6 
U 
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Exact 

Lyzwa's approach 
Our approach 



4 6 
U 



10 



-1 



t=0.3 




Exact 

Lyzwa's approach 
Our approach 



4 6 
U 
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10 



LU 



-1 



t=0.4 




Exact 

Lyzwa's approach 
Our approach 



4 6 
U 
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U=1, t f =0.2 



U=1, t f =0.4 




0.4 
0.35 
0.3 



4 0.25 



Exact 

Lyzwa's approach 
Our approach 




U=10, t=0.2 



U=10, t=0.4 




0.6 



0.4 



Exact 

- - Lyzwa's approach 
Our approach 




